//Converts pre-generated staggered prop to naive prop
//IF STAGGERED IN JIM FORMAT EXISTS, CAN USE BLANK FILE FOR ARGV[1], BLANK FOR ARGV[2], STAGGERED FILE FOR ARGV[3], AND NAME OF OUTPUT FILE FOR ARGV[4]
#include "/home/rtevans/fermiqcd/fermiqcd.h"
#include "/home/rtevans/fermiqcd/fermiqcd_JIM_IO.h"


//For staggered->naive propagator conversion
Matrix OmegaMatrix(site x) { 
  Matrix M(4,4); 
  M=1; 
  if(x(3)%2==1) M=Gamma[3]*M; 
  if(x(2)%2==1) M=Gamma[2]*M; 
  if(x(1)%2==1) M=Gamma[1]*M; 
  if(x(0)%2==1) M=Gamma[0]*M;  
  return M;
} 


int main(int argc, char** argv) {
  mpi. open_wormholes(argc,argv);
  //Staggered Prop to load
  char* info=argv[1]; 
  string configuration=argv[2];
  char gauge_filename [FILENAME_MAX];
  strncpy(gauge_filename,configuration.c_str(),FILENAME_MAX);

  string staggered_propagator=argv[3];
  char staggered_prop [FILENAME_MAX];
  strncpy(staggered_prop,staggered_propagator.c_str(),FILENAME_MAX);

  char* output=argv[4];
  //Set up lattice
  int nc=3;  
  int L[4];
  L[0]=(int) val(prompt(info,"T","48"));
  L[1]=(int) val(prompt(info,"X","16"));
  L[2]=(int) val(prompt(info,"Y","16"));
  L[3]=(int) val(prompt(info,"Z","16"));
  int nt=L[0];

  define_base_matrices("FERMILAB");

  generic_lattice mylattice(4,L,default_partitioning0,
			    torus_topology,0,3,false);
  mdp_matrix_field chi(mylattice,nc,nc);
  fermi_propagator S_l(mylattice,nc);
  site x(mylattice);
  
  //Create Omega Matrix for staggered to naive conversion
  mdp_matrix_field Omega(mylattice,4,4);
  forallsites(x) Omega(x)=OmegaMatrix(x); 
  
  coefficients staggered_quark;
  staggered_quark["mass"]=val(prompt(info,"mass",".0492"));
  staggered_quark["u0"]=val(prompt(info,"u0",".8564"));
  //////////Naive prop if staggered does not exist//////////////////
  if(staggered_propagator=="NONE"){
    
    mdp_real absolute_precision=val(prompt(info,"ABSOLUTE_PRECISION","1e-21"));
    mdp_real relative_precision=val(prompt(info,"RELATIVE_PRECISION","1e-6"));

    gauge_field U(mylattice,nc);
    gauge_field V(mylattice,nc);  

    //Use MILC config in JIM_IO format
    //read in gauge config, associate gauge config with U fields
    MetaHash UserMeta;
    JIM_load<float> ( U, gauge_filename, UserMeta );
    compute_em_field(U);
    
    //Improve staggered action, improved links are V
    mdp_real avg_plaq=average_plaquette(U);  
    lepage_improved_links(V,U,lepage_coefficients(pow(staggered_quark["u0"],4), "Full"));
    default_staggered_action=StaggeredAsqtadActionFast::mul_Q;  
    default_staggered_inverter=StaggeredBiCGUML::inverter; 
    staggered_field chi1(mylattice,nc); 
    staggered_field chi2(mylattice,nc); 
    //Naive quark
    for(int i=0;i<nc;i++) {
      chi1=0;
      //Set source
      if(on_which_process(mylattice,0,0,0,0)==ME) {
	x.set(0,0,0,0);
	chi1(x,i)=1;  
      }
      chi1.update();
      //Create Staggered prop
      mul_invQ(chi2,chi1,V,staggered_quark,1e-20,0,2000);
      //Convert to Naive prop
      forallsites(x) { 
	for(int j=0;j<nc;j++) 
	  for(int a=0;a<4;a++) 
	    for(int b=0;b<4;b++)
	      S_l(x,a,b,j,i)=Omega(x,a,b)*chi2(x,j);
      }
    }
    char outfile[1024];
    sprintf(outfile,"%s_%f_%s",output,staggered_quark["mass"],"naive");
    S_l.save(outfile); 
  }
  else{
    /////////////////////////////////////////////////////////////
    ///////////Staggered to naive if Staggered Exists///////////
    //Load Staggered Prop
    MetaHash StagMeta;
    JIM_load<float> (chi,staggered_prop, StagMeta);
    //Convert to Naive prop
    for(int i=0;i<nc;i++) {
      for(int j=0;j<nc;j++) 
	for(int a=0;a<4;a++) 
	  for(int b=0;b<4;b++)
	    forallsites(x) {
	    S_l(x,b,a,j,i)=Omega(x,b,a)*chi(x,i,j);
	  }
    }
    mdp<<"Existing Staggered has been Rotated" <<endl;
  char outfile[1024];
  sprintf(outfile,"%s_%f_%s",output,staggered_quark["mass"],"naive");
  S_l.save(outfile);
  }
  mpi.close_wormholes();
}


